#!/usr/bin/Rscript

time.offset.default = 0
mean.radius = 6371.009 # km
wind.dir.max.p = 0.01

gps2coords <- function(gps.df,
                       longitude.name, latitude.name, time.name,
                       time.offset = time.offset.default) {
    if (missing(gps.df))
        stop("missing argument: gps.df")
    if (missing(longitude.name))
        stop("missing argument: longitude.name")
    if (missing(latitude.name))
        stop("missing argument: latitude.name")
    if (missing(time.name))
        stop("missing argument: time.name")
    
    hms <- sub("^.*?(\\d{,2}:\\d{2}:\\d{2}(\\.\\d)?(\\s*[ap]m)?).*$", "\\1",
               gps.df[[time.name]], ignore.case=TRUE)
    is.pm <- grepl("\\s*pm$", hms, ignore.case=TRUE)
    hms <- sub("^(.*)\\s*[ap]m$", "\\1", hms, ignore.case=TRUE)
    time <- sapply(strsplit(as.character(hms), ":"), function(x) {
        x <- as.numeric(x); (x[1] * 60 + x[2]) * 60 + x[3]
    })
    time <- time + ifelse(is.pm, 12*60*60, 0)

    longitude <- gps.df[[longitude.name]]
    latitude <- gps.df[[latitude.name]]
    long.lat <- data.frame(time.sec = time,
                           long.rad = longitude * pi / 180,
                           lat.rad = latitude * pi / 180)

    temp <- data.frame(cos.lat = cos(long.lat$lat.rad),
                       pz=sin(long.lat$lat.rad))
    temp$px <- temp$cos.lat * cos(long.lat$long.rad)
    temp$py <- temp$cos.lat * sin(long.lat$long.rad)

    centroid <- as.data.frame(t(colMeans(temp)))
    mm <- with(centroid, px*px + py*py + pz*pz)

    centroid <- centroid / sqrt(mm)
    centroid$sin.lat <- centroid$pz
    centroid$cos.lat <- sqrt(1-centroid$pz^2)

    temp$p.c <- temp$px * centroid$px + temp$py * centroid$py + temp$pz * centroid$pz
    temp$rx <- (temp$py * centroid$px - temp$px * centroid$py) / temp$p.c / centroid$cos.lat
    temp$ry <- (temp$pz / temp$p.c - centroid$sin.lat) / centroid$cos.lat
    
    pos <- data.frame(time = long.lat$time.sec/3600 + time.offset,
                      x = temp$rx * mean.radius,
                      y = temp$ry * mean.radius)
    centroid.rad <- c(atan2(centroid$py, centroid$px),
                                   asin(centroid$sin.lat))
    names(centroid.rad) <- c("longitude", "latitude")
    attr(pos, "centroid.rad") <-  centroid.rad
    attr(pos, "longitude") <- longitude
    attr(pos, "latitude") <- latitude
    attr(pos, "time.offset") <- time.offset
    
    return(pos)
}

gps2coords.n82 <- function(df.n82, time.offset = time.offset.default) {
    gps2coords(df.n82,
               longitude.name = "longitude",
               latitude.name = "latitude",
               time.name = "real.time",
               time.offset = time.offset)
}

gps2coords.gpsprune <- function(df.gpsprune,
                                time.offset = time.offset.default) {
    gps2coords(df.gpsprune,
               longitude.name = "Longitude",
               latitude.name = "Latitude",
               time.name = "Time",
               time.offset = time.offset)
}


rowdelta <- function(df) {
    nr <- nrow(df)
    return(df[2:nr,] - df[1:(nr-1),])
}

rowrate <- function(df, delta = rowdelta(df)) {
    rate <- cbind(df[2:nrow(df), "time", drop=FALSE],
                  delta[, c("x","y")] / delta[, "time"])
    na.rows <- which(is.infinite(rate$x) | is.infinite(rate$y))
    if (length(na.rows) > 0)
        rate[na.rows, c("x", "y")] <- NA
    return(rbind(data.frame(time=df[1,"time"], x=NA, y=NA),
                 rate))
}

tan.predom.dir <- function(velocity) {
    model <- lm(y~x, velocity)
    summ <- summary(model)
    tan <- summ$coefficients["x", "Estimate"]
    attr(tan, "p.val") <- summ$coefficients["x", "Pr(>|t|)"]
    return (tan)
}

add.polar <- function(df) {
    df$r <- sqrt(df$x^2 + df$y^2)
    df$a <- atan2(df$y, df$x)
    return(df)
}

pos.plot <- function(pos,
                     vel = rowrate(pos),
                     wind.dir = calculate.wind.dir(vel),
                     planing = calculate.planing(pos,
                         velocity = vel,
                         wind.dir=wind.dir,
                         ...),
                     vlim = c(0,30), # km/h
                     colour.low = if (dark.foreground) hcl(120, 30, 80)
                                  else hcl(-30, 50,10),
                     colour.high = if (dark.foreground) hcl(0, 100, 10)
                                  else hcl(90, 90, 90),
                     show.map = NULL, # or "openstreetmap" ("OSM"), "google"
                     map.google.type = "hybrid",
                     dark.foreground = is.null(show.map) ||
                       show.map=="OSM" ||
                       map.google.type=="roadmap" || map.google.type=="terrain",
                     ...) {
    if (is.null(vel$r)) {
        vel <- add.polar(vel)
    }
    data <- data.frame(x=pos$x, y=pos$y,
                       planing = planing,
                       speed = vel$r)

    library(ggplot2)
    library(scales)
    
    plot <- ggplot(data) +
        aes(x=x, y=y, colour=speed, colour=speed, shape=planing, size=planing) +
        coord_fixed() +
        scale_x_continuous("West <--> East") +
        scale_y_continuous("South <--> North") +
        scale_shape_manual(values=c(8, 5, 1), drop=FALSE) +
        scale_size_manual(values=c(2, 2, 1), drop=FALSE) +
        scale_colour_gradient("speed (km/h)",
                              low=colour.low, high=colour.high,
                              lim=vlim, oob=squish)

    if (!is.null(show.map)) {
        longitude <- attr(pos, "longitude")
        latitude <- attr(pos, "latitude")
        centroid <- attr(pos, "centroid.rad") * 180/pi
        if (length(centroid) != 2)
            centroid <- c(longitude = mean(longitude),
                          latitude = mean(latitude))

        bg.map <- fetch.bg.map(longitude.range = extend.range(c(min(longitude),
                                   max(longitude))),
                               latitude.range = extend.range(c(min(latitude),
                                   max(latitude))),
                               centroid = centroid,
                               plot.x.range = extend.range(c(min(pos$x),
                                   max(pos$x))),
                               plot.y.range = extend.range(c(min(pos$y),
                                   max(pos$y))),
                               map.type = show.map,
                               ...)

        if (inherits(bg.map, "proto"))
            plot <- plot + bg.map
    }

    plot +
        geom_abline(slope = tan(wind.dir),
                    colour=if (dark.foreground) muted("green")
                           else hcl(120, 40, 60)) +
        geom_point()
}

extend.range <- function(range, by.ratio = 1.1) {
    mid <- (range[1] + range[2])/2
    width <- (range[2] - range[1]) * by.ratio
    c(mid-width/2, mid+width/2)
}

fetch.bg.map <- function(longitude.range, latitude.range,
                         centroid,
                         plot.x.range, plot.y.range,
                         map.type,
                         map.google.zoom,
                         map.google.scale = 2,
                         map.google.type = "hybrid", # or "terrain", "satellite", "roadmap"
                         map.osm.scale.per.km = 8000,
                         ...) {
    library(ggmap)
    if (map.type == "openstreatmap" || map.type == "OSM") {
        map.bbox <- c(left=longitude.range[1], bottom=latitude.range[1],
                      right=longitude.range[2], top=latitude.range[2])
        map.scale <- (plot.x.range[2]-plot.x.range[1]) *
            map.osm.scale.per.km

        map <- get_openstreetmap(bbox = map.bbox, scale=round(map.scale),
                                 messaging = TRUE)

        return(annotation_raster(map,
                                 xmin=plot.x.range[1],
                                 xmax=plot.x.range[2],
                                 ymin=plot.y.range[1],
                                 ymax=plot.y.range[2]))
    } else if (map.type == "google") {
        long.width <- longitude.range[2] - longitude.range[1]
        lat.width <- latitude.range[2] - latitude.range[1]
        x.width <- plot.x.range[2] - plot.x.range[1]
        y.width <- plot.y.range[2] - plot.y.range[1]

        if (missing(map.google.zoom))
            zoom <- 15 + log2(3) - log2(1.1*max(x.width, y.width))
        else
            zoom <- map.google.zoom
        zoom <- max(0, min(floor(zoom), 21))
        map <- get_googlemap(center = centroid,
                             zoom=zoom, scale=map.google.scale,
                             maptype=map.google.type)

        bb <- attr(map, "bb")
        xmin <- (bb$ll.lon - centroid["longitude"])/long.width * x.width
        xmax <- (bb$ur.lon - centroid["longitude"])/long.width * x.width
        ymin <- (bb$ll.lat - centroid["latitude"])/lat.width * y.width
        ymax <- (bb$ur.lat - centroid["latitude"])/lat.width * y.width

        return(annotation_raster(map,
                                 xmin=xmin, xmax=xmax,
                                 ymin=ymin, ymax=ymax))
    }
    return(NA)
}

pos.anim <- function(pos,
                     vel = rowrate(pos),
                     wind.dir = calculate.wind.dir(vel),
                     planing = calculate.planing(pos,
                         velocity = vel,
                         wind.dir=wind.dir,
                         ...),
                     output.filename.prefix,
                     show.map = NULL, # or "openstreetmap" ("OSM"), "google"
                     map.google.type = "hybrid",
                     dark.foreground = is.null(show.map) ||
                         show.map=="OSM" ||
                         map.google.type=="roadmap" ||
                         map.google.type=="terrain",
                     ...) {
    if (is.null(vel$r)) {
        vel <- add.polar(vel)
    }
    data <- data.frame(x=pos$x, y=pos$y, planing = planing,
                       speed = vel$r, time = pos$time)

    library(ggplot2)
    library(scales)
    
    longitude <- attr(pos, "longitude")
    latitude <- attr(pos, "latitude")
    centroid <- attr(pos, "centroid.rad") * 180/pi
    if (length(centroid) != 2)
        centroid <- c(longitude = mean(longitude),
                      latitude = mean(latitude))

    longitude.range = extend.range(c(min(longitude), max(longitude)))
    latitude.range = extend.range(c(min(latitude), max(latitude)))
    plot.x.range = extend.range(c(min(pos$x), max(pos$x)))
    plot.y.range = extend.range(c(min(pos$y), max(pos$y)))
    plot.speed.range = extend.range(c(min(vel$r, na.rm=TRUE),
        max(vel$r, na.rm=TRUE)))

    point.colours = c(
        no        = hcl(h=60,  c=80,   if (dark.foreground) 40 else 60),
        starboard = hcl(h=120, c=90, l=if (dark.foreground) 20 else 80),
        port      = hcl(h=0,   c=90, l=if (dark.foreground) 20 else 80))

    plot <- ggplot(data.frame()) +
        aes(x=x, y=y, colour=planing, size=speed) +
        coord_fixed() +
        scale_x_continuous("West <--> East", limits = plot.x.range) +
        scale_y_continuous("South <--> North", limits = plot.y.range) +
        scale_colour_manual(values=point.colours, guide = FALSE) +
        scale_size_area(limits=plot.speed.range, guide=FALSE)

    if (!is.null(show.map)) {
        bg.map <- fetch.bg.map(longitude.range = longitude.range,
                               latitude.range = latitude.range,
                               centroid = centroid,
                               plot.x.range = plot.x.range,
                               plot.y.range = plot.y.range,
                               map.type = show.map,
                               ...)

        if (inherits(bg.map, "proto"))
            plot <- plot + bg.map
    }

    plot <- plot +
        geom_abline(slope = tan(wind.dir),
                    colour=if (dark.foreground) muted("green")
                           else hcl(120, 40, 60)) +
        geom_point()

    png(paste0(output.filename.prefix, "-%04d.png"), 800, 600)
    for (i in 1:nrow(data)) {
        if (i %% 30 == 0)
            message("generating frame ", i)
        d <- data[i, ]
        p <- plot %+% d +
            ggtitle(paste(time2hms(d$time),
                          paste0(if (is.na(d$speed)) "    "
                                 else sprintf("%4.1f", d$speed),
                                 "km/h"),
                          sep = "\t"))
        print(p)
    }
    dev.off()
}


calculate.planing <- function(pos,
                              vel.cross = calculate.vel.cross(pos, ...),
                              planing.speed.min = 15,
                              ...) {
    factor(ifelse(vel.cross <= -planing.speed.min, "starboard",
                  ifelse(vel.cross >= +planing.speed.min, "port",
                         "no")),
           levels=c("starboard", "port", "no"))
}

calculate.wind.dir <- function(vel) {
    cotan.wind <- tan.predom.dir(vel)
    wind.dir <- atan(cotan.wind) + pi/2
    if (is.null(vel$r)) {
        vel = add.polar(vel)
    }
    cosines <- (vel$x * cos(wind.dir) + vel$y * sin(wind.dir)) / vel$r
    if (sum(cosines, na.rm=TRUE) > 0)
        return(wind.dir + pi)
    return(wind.dir)
}

calculate.vel.cross <- function(pos,
                                velocity = rowrate(pos),
                                wind.dir = calculate.wind.dir(vel)) {
    velocity$x * sin(wind.dir) - velocity$y * cos(wind.dir)
}

cluster <- function(pos,
                    planing = calculate.planing(pos,
                        planing.speed.min = planing.speed.min, ...),
                    planing.speed.min = 15, #km/h
                    gap.time.max = 5/3600, # hours
                    gap.vel.min = planing.speed.min*0.9, #km/h
                    ...) {
    n = length(planing)
    q = integer(length=n)
    if (n < 1) return(q)
    gap.after <- NA
    gap.i <- NA
    cluster <- "no"
    cluster.id <- 0
    for (i in 1:n) {
        curr <- planing[i]
        if (is.na(curr)) curr <- "no"
        if (curr == cluster) {
            if (cluster == "no") {
                q[i] <- NA
                next
            } else {
                q[i] <- cluster.id
                next
            }
        } else {
            if (curr == "no") {
                gap.after <- cluster
                gap.i <- i
                cluster <- "no"
                q[i] <- NA
                next
            } else {
                if (!is.na(gap.after)
                    && gap.after == curr
                    && pos[i, "time"] - pos[gap.i-1, "time"] <= gap.time.max
                    && sqrt( (pos[i, "x"] - pos[gap.i-1, "x"])^2
                            +(pos[i, "y"] - pos[gap.i-1, "y"])^2) /
                          (pos[i, "time"] - pos[gap.i-1, "time"]) >=
                       gap.vel.min) {
                    # join gap
                    for (j in gap.i:(i-1)) q[j] = cluster.id
                    gap.after <- NA
                    gap.i <- NA
                    cluster <- curr
                    q[i] <- cluster.id
                    next
                } else {
                    gap.after <- NA
                    gap.i <- NA
                    cluster <- curr
                    cluster.id <- cluster.id + 1
                    q[i] <- cluster.id
                    next
                }
            }
        }
    }
    return(q)
}

analyse.clusters <- function(pos,
                             cl=cluster(pos, ...),
                             time.interval.min = 5/3600, # hr
                             vel.drop.max = 0.95,
                             ...) {
    output <- data.frame(time=numeric(0),
                         duration=numeric(0),
                         speed=numeric(0),
                         cluster=numeric(0))
    for (c in levels(factor(cl))) {
        pvc1 <- pos[!is.na(cl) & cl==c, , drop=FALSE]
        n <- nrow(pvc1)
        if (n < 2) next
        v.max <- -Inf
        interval <- NULL
        for (i in 1:(n-1)) {
            diff.time <- pvc1$time[(i+1):n] - pvc1$time[i]
            diff.v <- sqrt((pvc1$x[(i+1):n] - pvc1$x[i])^2 +
                           (pvc1$y[(i+1):n] - pvc1$y[i])^2) / diff.time
            m <- length(diff.time)
            j <- 1
                while (j <= m && diff.time[j] < time.interval.min)
                    j <- j+1
                if (j > m) break;
                v = diff.v[j]
                j <- j+1
                while (j <= m && diff.v[j] >= v * vel.drop.max) {
                    j <- j+1
                }
                j0 <- j-1
                v = diff.v[j0]
                if (v > v.max) {
                    v.max <- v
                    interval <- rownames(pvc1)[c(i, i+j0)]
                }
        }
        if (!is.null(interval)) {
            p <- pos[interval,]
            dur <- p[2, "time"] - p[1, "time"]
            output.row <- data.frame(time=p[1, "time", drop=FALSE],
                                     duration=dur,
                                     speed=v.max,
                                     cluster=c)
            output <- rbind(output, output.row)
        }
    }
    output
}

time2hms <- function(t) { # t in hr
    h <- as.integer(t)
    m0 <- (t-h) * 60
    m <- as.integer(m0)
    s0 <- (m0-m) * 60
    s <- as.integer(round(s0))
    sprintf("%02d:%02d:%02d", h, m, s)
}


top.speeds.plot <- function(pos = NULL,
                            cl.an = analyse.clusters(pos, ...),
                            ...) {
    an <- cl.an
    if (nrow(an)<1) return(NA)
    an$time <- an$time

    library(ggplot2)
    library(scales)

    ggplot(an) +
        aes(x=duration*3600, y=speed, label=cluster, colour=time) +
        geom_point() +
        geom_text(show_guide=F, vjust=-0.4) +
        scale_x_continuous("duration (s)") +
        scale_y_continuous("speed (km/h)") +
        scale_colour_gradientn(name="time (h:m:s)", labels=time2hms,
                               colours=hcl(seq(0,360,by=60),100,80)) +
        theme(panel.background=element_rect(fill="black"),
              panel.grid.major=element_line(colour=grey(.5)),
              panel.grid.minor=element_line(colour=grey(.3)),
              legend.background=element_rect(fill="black"),
              legend.title=element_text(colour="white"),
              legend.text=element_text(colour="white")
              )
}


top.speeds2.plot <- function(pos = NULL,
                             cl.an = analyse.clusters(pos, ...),
                             ...) {
    an <- cl.an
    if (nrow(an)<1) return(NA)

    today <- as.POSIXlt(Sys.time(), tz="GMT")
    today$hour <- 0
    today$min <- 0
    today$sec <- 0
    an$time.lt <- today + an$time * 3600

    library(ggplot2)
    library(scales)

    ggplot(an) +
        aes(colour=duration*3600, y=speed, label=cluster, x=time.lt) +
        geom_point() +
        geom_text(show_guide=F, vjust=-0.4) +
        scale_x_datetime("time") +
        scale_y_continuous("speed (km/h)") +
        scale_colour_continuous(name="duration (s)",
                                low=gray(0.3),
                                high=col2hcl("yellow", c=100, l=90)) +
        theme(panel.background=element_rect(fill="black"),
              panel.grid.major=element_line(colour=grey(.5)),
              panel.grid.minor=element_line(colour=grey(.3)),
              legend.background=element_rect(fill="black"),
              legend.title=element_text(colour="white"),
              legend.text=element_text(colour="white")
              )
}


analyse <- function(filename, # either one
                    format.subtype = "n82", # or "gpsprune"
                    df, # pre-fetched data frame
                    time.offset = time.offset.default,
                    wind.dir, # for overriding the computed wind dir.
                    ...) {
    if (missing(df)) {
        if (missing(filename))
            stop("neither 'df' nor 'filename' is given")
        if (grepl("\\.csv", filename))
            df <- read.csv(filename)
        else
            df <- read.delim(filename)
    } else {
        if (!missing(filename))
            warning("'filename' ignored since 'df' is given")
    }

    if (format.subtype == "n82")
        pos <- gps2coords.n82(df, time.offset = time.offset)
    else if (format.subtype == "gpsprune")
        pos <- gps2coords.gpsprune(df, time.offset = time.offset)
    else if (format.subtype == "gpx2txt")
        pos <- gps2coords.gpsprune(df, time.offset = 0)
    else
        stop("invalid 'format.subtype': ", format.subtype,
             "; acceptable: 'n82' or 'gpsprune'")

    pos.delta <- rowdelta(pos)
    vel <- add.polar(rowrate(pos, pos.delta))

    if (missing(wind.dir)) {
        wind.dir <- calculate.wind.dir(vel)
        if (attr(wind.dir, "p.val") > wind.dir.max.p)
            warning("big p-value for wind direction: ",
                    attr(wind.dir, "p.val"))
    }

    vel.cross = calculate.vel.cross(pos, vel, wind.dir)
    planing = calculate.planing(vel.cross=vel.cross, ...)
    turns <- count.turns(planing)

    cl <- cluster(pos=pos, planing=planing, ...)

    an <- analyse.clusters(pos, cl=cl, ...)
    an$hms <- time2hms(an$time)
    an$dur <- an$duration * 3600

    time.range <- range(pos$time, na.rm=T)
    distance <- sum(sqrt(pos.delta$x^2 + pos.delta$y^2))
    
    results = NULL

    results$pos <- pos
    results$duration <- time.range[2] - time.range[1]
    results$distance <- distance
    results$velocity <- vel
    results$wind.dir <- wind.dir
    results$planing <- planing
    results$turns <- turns
    results$cl <- cl
    results$cl.an <- an
    results$time.offset <- time.offset

    attr(results, "class") <- "analysis"
    return(results)
}

count.turns <- function(planing) {
    planing[planing == "no"] <- NA
    planing <- na.omit(planing)
    n <- length(planing)
    if (n < 2) return(0)
    diff <- planing[1:(n-1)] != planing[2:n]
    return(length(which(diff)))
}

print.analysis <- function(analysis) {
    wind.bearing <- 90 - analysis$wind.dir * 180 / pi
    while (wind.bearing < 0) {
        wind.bearing <- wind.bearing + 360
    }
    cat("points: ", nrow(analysis$pos), "\n", sep="")
    cat("duration: ", time2hms(analysis$duration), "\n", sep="")
    cat("distance: ", analysis$distance, " km\n", sep="")
    cat("centroid: ", attr(analysis$pos, "centroid.rad") * 180 / pi, "\n")
    cat("wind bearing: ", wind.bearing, "\n", sep="")
    cat("turns: ", analysis$turns, "\n", sep="")
    cat("planing count: ", max(c(0,analysis$cl), na.rm=TRUE), "\n", sep="")
    cat("planing clusters: ", nrow(analysis$cl.an), "\n", sep="")
    cat("\nPlanings by speed\n")
    print(head(analysis$cl.an[order(analysis$cl.an$speed, decreasing=TRUE), ],
               10))
    cat("\nPlanings by duration\n")
    print(head(analysis$cl.an[order(analysis$cl.an$duration, decreasing=TRUE),],
               10))
    invisible(analysis)
}

plot.analysis <- function(analysis, ...) {
    longlat <- data.frame(longitude = attr(analysis$pos, "longitude"),
                          latitude  = attr(analysis$pos, "latitude"))
    invisible(list(pos = pos.plot(analysis$pos,
                       velocity = analysis$velocity,
                       wind.dir = analysis$wind.dir,
                       planing = analysis$planing,
                       ...),
                   top.speeds = top.speeds.plot(cl.an = analysis$cl.an,
                       ...),
                   top.speeds2 = top.speeds2.plot(cl.an = analysis$cl.an,
                       ...)
                   ))
}


cmd.analyse <- function(f, m, turn.around=FALSE) {
if (!is.na(f)) {
    r <- tryCatch(analyse(f),
             error = function(e) { analyse(f, "gpx2txt") })
    if (turn.around)
        r <- tryCatch(analyse(f, wind.dir = r$wind.dir+pi),
                      error = function(e) {
                          analyse(f, "gpx2txt", wind.dir=r$wind.dir+pi)
                      })
    print(r)
    if (is.na(m))
      p <- plot(r)
    else
      p <- plot(r, show.map=m)

    f.basename <- sub("[.][^.]*$", "", basename(f))
    out.pdf.name <- paste0(f.basename, ".pdf")
    out.txt.name <- paste0(f.basename, ".analysis.txt")

    pdf(out.pdf.name)
    print(p$pos)
    if (!is.na(p$top.speeds))
      print(p$top.speeds)
    if (!is.na(p$top.speeds2))
      print(p$top.speeds2)
    z <- dev.off()
    message("plotted: ", out.pdf.name)

    message("saving analysis report to: ", out.txt.name)
    sink(out.txt.name)
    print(r)
}
}

cmd.anim <- function(f, m, o) {
    r <- tryCatch(analyse(f),
                  error = function(e) { analyse(f, "gpsprune", time.offset=8) })

    pos.anim(r$pos, show.map=m, output.filename.prefix=o)
}

usage <- function() {
    message("Sub-commands:")
    message("\tanalyse filename google|OSM")
    message("\tanim filename google|OSM output_filename_prefix")
}

args <- commandArgs(trailingOnly=TRUE)
c = args[1]
if (is.na(c)) {
    usage()
} else {
    if (c == "analyse") {
        cmd.analyse(args[2], args[3])
    } else if (c == "analyse180") {
        cmd.analyse(args[2], args[3], turn.around=TRUE)
    } else if (c == "anim") {
        cmd.anim(args[2], args[3], args[4])
    } else {
        usage()
        stop("unknown command '", c, "'; valid commands are: analyse anim")
    }
}

#end
